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We investigate the effects of strong number fluctuations on traveling waves in the Fisher- 
Kolmogorov reaction-diffusion system. Our findings are in stark contrast to the commonly used 
deterministic and weak-noise approximations. We compute the wave velocity in one and two spatial 
dimensions, for which we find a linear and a square-root dependence of the speed on the particle 
density. Instead of smooth sigmoidal wave profiles, we observe fronts composed of a few rugged kinks 
that diffuse, annihilate, and rarely branch; this dynamics leads to power-law tails in the distribution 
of the front sizes. 
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Traveling waves are a common phenomenon in many 
systems that combine diffusion and reaction of particles. 
Familiar examples are the combustion fronts running 
through a premixed reactive gas, or the expanding fronts 
of bacterial colonies. The standard model for these wave- 
like phenomena is the stochastic Fisher-Kolmogorov- 
Petrovsky-Piscounov (sFKPP) equation, the simplest 
non-linear model that blends diffusion, growth and num- 
ber fluctuations. It has been used widely in population 
genetics [l[, ecology 0, epidemiology [5], chemical ki- 
netics Q, and recently even in quantum chromodynam- 
ics 0. 

Our current understanding of traveling waves in reac- 
tion diffusion systems is shaped by studies of the sFKPP 
equation that have focused on the weak-noise regime. 
The opposite limit of strong noise is relatively unex- 
plored, yet arguably of equal importance. Not only does 
it occur naturally, when the driving forces of traveling 
waves are weak [5[; it also unravels fundamentally dif- 
ferent and often counterintuitive aspects of the sFKPP 
equation. The focus of this letter is on two intimately 
related questions of how fast these traveling waves move 
and what their shape is. 

Recently, the velocity of a one-dimensional wave has 
been computed in the strong noise limit with the help 
of a duality between the sFKPP equation and A ^ 
A + A reaction-diffusion system [Q]. Here, we treat the 
growth (reaction) term as a small perturbation, and con- 
struct a more direct and more powerful method to calcu- 
late the speed of the wave than the one used in Ref . Q . 
Our technique has three important advantages. First, 
the perturbation expansion enables us to compute the 
wave velocity in two dimensions. We find a square root 
dependence of the terminal velocity on the particle den- 
sity. Second, our approach allows studying the shape of 



the front on time and length scales that are dominated by 
the noise. We find that, in one dimension, the noise alone 
can stabilize the front and limit its diffusive broadening. 
The front size distribution exhibits power-law tails due 
to spontaneous creation and subsequent annihilation of 
new transition regions. Third, the perturbation expan- 
sion can be applied to more general models with polyno- 
mial reaction terms of higher order. 

Without any loss of generality, we discuss the sFKPP 
equation from the point of view of population genetics, 
where it is easy to interpret, simulate, and potentially 
test [Ij. In this context, the sFKPP is used to describe 
how a mutation that increases the growth rate of its car- 
rier spreads through a homogeneous population [H (see 
also Ref. 0] for a recent review). In one dimension, the 
relative abundance p(t, x) S [0, 1] of a beneficial mutation 
at time t and position x is described by 

% ^ + - P) + - P)ri{t, x). (1) 

The diffusion term is due to the short-range migration 
of the individuals carrying the mutations. The diffusiv- 
ity D is proportional to the average dispersal distance in 
one generation. The reaction term, ap{l — p), accounts 
for different fitnesses of the mutant and the "wild" type, 
and the reaction rate a > is the difference in their 
growth rates. The last term on the right hand side of 
Eq. ^ describes the sampling error during reproduc- 
tion, and is commonly referred to as genetic drift or 
number fluctuations. The strength of the noise 6 > 
is inversely proportional to the population density, and 
the (Ito) white noise r]{t, x) satisfies the following condi- 
tion: {r]{ti,xi)ri{t2,X2)) = 5{ti-t2)5{xi- X2), vfYieve 5{-) 
stands for Dirac's delta function. 

One of the most important predictions of the sFKPP 
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equation is the velocity v of an isolated wave that moves 
from the left half-space occupied by mutants into the 
right half-space occupied by the less fit wild type. The 
boundary between the half-spaces is assumed to be sharp 
initially, p{0,x) = 1 — 0{x), where 6{x) is the Heaviside 
step function. In the deterministic limit (6 = 0), the wave 
acquires a stationary shape with exponential tales, and 
the velocity of the front approaches vp = 2\jDa [H, 
Surprisingly, even weak noise, &^ <C aZ), gives rise to large 
corrections to the velocity, v = vp — 0[ln~^(l/6)] Q, in- 
dicating the importance of the noise. Here, we access the 
strong noise regime, \? ^ aD, by constructing a pertur- 
bation expansion in the reaction rate a, while treating 
the noise term exactly. 

Let us define the instantaneous velocity of the wave 
as the average growth rate of the mutants, v = 

(^J^^p(t,x)dx'j, which is consistent with the usual 

definition in the deterministic limit. We then take the 
time derivative inside the integral and use Eq. ([T]) to 
eliminate after integrating by parts and averaging, 
we obtain an alternative expression for the velocity, 



— CL (p{t, x)[l ~ p{t, x)]) dx, 
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which relates the wave speed and the instantaneous wave 
profile encoded in the dynamical field p{t, x). 

Note that, to the first order in a, the instantaneous 
velocity is given by a/(t)/2, where the moment I{t) = 
/_oo (2p(^i — p{t, x)]) dx\a=o is evaluated in the neu- 
tral limit a — 0, when neither the mutant nor the wild 
type have a selective advantage. The neutral model is 
exactly solvable because the hierarchy of the moment 
equations closes 0,[3- For the purpose of this paper, it is 
sufficient to consider only the two-point correlation func- 
tion H(t, Xi,X2) = {pit, Xi)[l - p{t, X2)]) + {p{t, X2)[l - 

p{t,Xi)]), which is known in population genetics as the 
average spatial heterozygosity. H{t, xi,X2) is the aver- 
age probability that, at a time t, two individuals sam- 
pled at xi and X2 carry different genetic variants. The 
equation of motion for H is obtained by differentiating its 
definition with respect to time and eliminating ^ with 
the help of Eq. ([T|) . Note that the rules of the Ito calculus 
must be used to properly account for the effects of the 
noise 0, . The result is 
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Equation ^ can be understood intuitively as follows. 
The probability of sampling two different genetic vari- 
ants can be traced to the initial condition by following 
the lineages of the sampled variants backward in time. 
Then, H(t, xi, X2) changes due to the diffusion and due 
to the coalescence of the lineages represented by the term 



proportional to the delta function. The lineages coalesce 
if the genetic variants have a common ancestor, which can 
happen only when the lineages occupy the same point in 
space at the same time. 

We compute H(t, xi,X2) because it has information 
about the shapes of the wave front and is related to I{t), 
and thus v, by I{t) — H{t,x,x)dx. To this end, 
we introduce new spatial variables, ^ = [xi + 0:2 )/2 
and X = xi — a;2, encoding the average position of two 
sampling points and the distance between them, respec- 
tively. After a Laplace transform in t and Fourier trans- 
form in ^, Eq. ([3]) can then be solved for H{s,k,x) = 
dte-"^ die-^^^H{t, ^, x). The solution reads 
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From the definition of /(<), it follows that the Laplace 
transform of I{t) equals i/(s, 0,0). Performing the in- 
verse Laplace transform, we get Vmvt^ao I {i) = iDb'^^. 
Equation ([2]) then implies that the terminal velocity 
equals 2aD/b in agreement with Ref. Q. 

Since only the regions with p{t, x) 0, 1 contribute 
to /(t), the finite limit of I{t) as t ^ 00 may look 
counterintuitive because it suggests a finite length of the 
transition regions in the neutral front; in other words, 
the noise term alone limits the diffusive widening of the 
front. Qualitatively, this phenomenon can be understood 
by noticing that the probability of reaching local fixa- 
tion (p = or p = I) is very large at the tails of the front, 
so the stationary shape of the front might be maintained 
by a balance of the diffusive spreading of the genetic vari- 
ants into each other's territory and their subsequent loss 
due to number fluctuations (genetic drift). 

We analyze this peculiar phenomenon further via parti- 
cle simulations. Since our perturbation analysis suggests 
that, on sufficiently small length scales, the wave dynam- 
ics can be well approximated by neglecting Darwinian 
selection, we set a = in the simulations. A snapshot of 
a typical transition region between p = I and p = is 
shown in the inset of Fig. [TJ the transition occurs on a 
very short length scale set by D/b. The inset also shows 
the local heterozygosity /io(i, C) = 2p(t, C)[l — p{t,0]i 
which is nonzero only at the kink. The new coordinate C 
is defined such that J^^ho{tX)d( = ho{tX)dC, i.e. 
point C = is always in the center of the wave. Such 
a definition allows us to focus on the shape of the wave 
by eliminating the diffusion of the front. It is then in- 
structive to characterize the average shape of the front 
by :P(C) = limt^oo(p(i, C)) and /C(C) = limt^oo(^o(i, C)>- 
In contrast to the narrow boundary shown in the inset 
of Fig. [U these average characteristics show power-law 
tails with exponents close to — f and —2 respectively; see 

Fig.m 
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FIG. 1: (Color online) Simulation of the neutral one- 
dimensional stepping stone model [lol ] . The model considers a 
line of sites (demes) labeled by an integer j. Each of the demes 
is occupied by A'^ individuals that carry one of the two genetic 
variants. Every generation, neighboring demes exchange Nm 
migrating organisms. Migration is followed by Wright-Fisher 
reproduction, i.e. the next generation is formed by A'^ organ- 
isms sampled from the Bernoulli distribution with the proba- 
bility of sampling a particular genetic variant equal to its frac- 
tion in the current generation. Here, — 10, and m = 0.1; 
the total number of demes is 1501. We ensure that the wave 
front remains within the simulated habitat by moving C = 
to the center of the habitat every third generation. The inset 
shows a snapshot of the wave front p{t, Q [red (darker) color] 
and the local heterozygosity /io(i,C) (solid line) after lO'^ gen- 
erations starting from the step function initial condition. The 
main plot shows J-{(^) (dashed curve) and IC{Q (solid curve) 
obtained over 10* generations. For large ^, the functions have 
the following asymptotic scalings J-{Q oc and A3(C) oc 
shown by dotted-dashed and dotted lines respectively. 



The power-law tails of /C(C) and J-{C) can also be 
inferred from the exact solution for the two-point cor- 
relation function given by Eq. One can see 
that H{s,k,x) = G{s,k)£{s,k,x), where G = l/{s + 
Dk^) is the diffusion Green's function, which describes 
the motion of the center of the wave, and £ describes 
the evolution of the shape of the wave front. This 
decomposition implies that the front diffuses with dif- 
fusivity D independent of the noise strength b set by 
the population density. Furthermore, since factoriza- 
tion in the Fourier and Laplace domains corresponds 
to convolution in the space and time domains, we can 
think of S{t',^',x) as a contribution to H{t,£^,x) from 
a wave that was present a distance £, — S,' away and 
time t — t' ago. Power laws in the front size distri- 
bution should be reflected in the asymptotic behavior 



Eq. (HD, where $(t) = j°°^d(,£{t,^,x = 0) - i'^/^ for 
large times. The new function $(t) is the contribution 
to H from a transition region present t ago anywhere 
in space. Since dt£{t,£,,x = 0) ~ for large ^, 
we may conclude that the probability that a kink is lo- 
cated ^ away from the center of the wave should decay 
like a power law with exponent —2. Indeed, this matches 
the observed behavior of /C(C) for large arguments. 

The algebraic decay of correlations inside a wave front 
and the slow approach of I{t) to its limit as t —* oo 
are in contrast to the narrowness of a typical wave front 
shown in the inset of Fig. [1] It is unlikely that a single 
narrow front relaxes so slowly, but several diffusing fronts 
can exhibit very slow relaxation, give rise to power-law 
tails of /C and T, and still have a finite value of I{t = 
oo). Indeed, we find spontaneous creation of new kinks 
in our simulations. This process can be interpreted as the 
following reaction: A {2z + 1)A, where A represents a 
kink, and z is an integer. Note that the reaction with z = 
1 is the most frequent. Neighboring kinks can also merge 
and subsequently disappear, which is equivalent to the 
annihilation reaction ^ -I- ^ ^ 0. Thus the behavior of 
the front can be described by the dynamics of Branching 
Annihilating Random Walks (BARW) 

Earlier studies of BARW with an even number of off- 
spring jllj l have found that, in one dimension, the parti- 
cles do not proliferate regardless of the birth rate, which 
is consistent with the finite value of I{t = oo). More- 
over, we can understand the asymptotic behavior of /C(C) 
and J-{C) by considering the dynamics of only three an- 
nihilating random walks ( ARW) ; higher number of ARW 
lead to subleading corrections as one can easily show by 
generalizing our analysis of three ARW. 

The statistical properties of ARW have been reviewed 
by Fisher [l^ . Here, we extend his analysis to compute 
the average number of ARW present at position which 
is proportional to /C(C) for large C- Note, from the ARW 
interpretation, it follows that J-{C) oc / IC{()d( oc ^"^ 
because, for large (, J-{C) is proportional to the proba- 
bility that the farthest transition region is at least C away 
from the origin. So, we only have to calculate /C(C)- 

Since three ARW eventually annihilate, and the pro- 
cess repeats, it is sufficient to consider only one cycle 
from A — > 3A to 3A — * A. Let P{t,xi,X2,X3) be 
the probability to find three ARW at time t at posi- 
tions xi > X2 > X3, which obeys a three-dimensional 
diffusion equation, with the following absorbing bound- 
ary conditions P{t, xi,X2,X2) ~ P{t,xi,xi,X'i) — 
P{t, xi,X2,xi) = 0. Upon using the solution of this 



I4I, we obtain /C(C) 



of £(t,e,x)- We get £{t,tx = 0) 



diffusion problem from Ref. 

J dtdxidx2dx3S[C ~ {xi - ^2i±i^±^)]P(i,a;i,a;2,X3) oc 
C^^. Also note that the survival probability of 
three ARW / dxidx2dx3P{t,xi,X2,X3) ~ t^^^^ for long 
times 12l|. which matches the behavior of ^(t). Thus, on 
large length scales, the one-dimensional sFKPP equation 
from seems indeed equivalent to a one-dimensional reaction- 
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diffusion system of BARW. 

What happens in higher dimensions? The noisy FKPP 
equation Eq. ((!]) can easily be extended to two spatial co- 
ordinates upon substituting x ^ r — {x,y), provided one 
limits the short-wave length variations of the field r) 
by a cutoff I, e.g., representing the lattice constant in the 
model. The speed of a planar wave front traveling in the 
X direction again takes the form of Eq. ^ , and depends 
on the moment I{t). Solving the two-dimensional neutral 
version of Eq. ^ yields 



/(.) = 



1 



(2d) 



(5) 



for the Laplace transform of I{t). Up to logarithmic cor- 
rections, which are typical for two dimensional diffusion 
problems [l^, we thus find that the second moment in- 
creases as I{t) ~ D^^^b~^t^^'^ for long times. Conse- 
quently, the wave speed should increase without bound as 
v{t) = aI{t)/2 cx t^/^, rendering the perturbation expan- 
sion singular. We circumvent this difficulty by noticing 
that, when the front is moving, it does not have enough 
time to fully relax as predicted by the neutral dynamics. 
Relaxation only occurs over a limited time t* , which is 
roughly the time when deterministic motion of the front 
is of the same order as its diffusive motion: vt* ~ \/2Dt*, 
i.e. t* — 2Dv~^. If this time is large, we may impose 
a self-consistency condition, v = 2al{t*), which leads 
to the scaling v ~ D^fajh (up to log-corrections) for 
the wave velocity. These heuristic arguments, which 
can be formalized by a multiple-scale ansatz [14|, thus 
predict a crossover from the no- noise speed v ^ 2\faD 
to V ^ D^fajh for large values of the noise strength, 
6 3> -D. Simulated (genetic) wave fronts indeed exhibit 
this crossover, see Fig. [21 

In conclusion, we have presented a perturbative ap- 
proach to the one-dimensional sFKPP reaction-diffusion 
system with strong number fluctuations, which can be 
extended to higher dimensions using self-consistency ar- 
guments. In one and two dimensions, we found a linear 
and a novel square-root relationship between the speed 
of traveling waves and the particle density, respectively. 
The wave profiles have also been analyzed in one dimen- 
sion using three different approaches: Simulations, the 
correspondence between BARW and sFKPP equation, 
and the exact solution for the two-point correlation func- 
tion. All of these approaches predict that the front size 
distribution has a power-law tail with a cutoff, which 
is in contrast to the exponential tails of deterministic 
Fisher waves. The power-law tail is best understood as a 
consequence of spontaneous creation of several transition 
regions that behave as branching annihilating random 
walks. Finally, we note that our perturbation expansion 
also applies to generalized sFKPP equations with higher 
order polynomial reaction terms. This can be used, for 




FIG. 2: (Color online) The speed of planar fronts in the two- 
dimensional stepping stone model. The simulations are iden- 
tical to the ones described in Fig.[l]except for two differences. 
First, we consider a (co- moving) square lattice of 250 x 250 
sites instead of a linear array. Second, a mutant is by a factor 
1 + ar more likely to reproduce than a wild type in a single 
generation time r. The simulations are initialized by a sharp 
step-hke front profile and run for 2 x 10^ generations. In the 
scaling regimes, the wave speeds closely follow a power law 
cx IT^^"^ — {Nni)^l'^ . Deviations are due to log-corrections, 
and mimic a power law with exponent slightly smaller than 
— 1/2. However, the scaling regimes almost collapse in a plot 
V = [v / V f) \/(iJrn versus h * m/a = {Ns)~^ , see the inset. 
The solid line is the solution of = c\* Ns/ log[c2 * u] with 
the fitting parameters ci — 0.45 and C2 — 1. 



instance, to show that the spread of recessive and domi- 
nant beneficial mutations is the same 
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